# This is the job submission script to resample the AWAP ASCII files of daily max and min temp
# To 250 m resolution for the WT plus Buffer regions

# Establish in directory

in.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPASCII/RawAWAP/'
setwd(in.dir)

# Load necessary libraries
necessary=c('SDMTools','sp','rgdal')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

	# Start by creating an object with row/col positions for the WT plus buffer regions
	# Then get XY's to obtain easting/northing
	# Then convert easting/northing to lat/long, storing in the same data frame
	# Then write out this file as a .csv to be imported for each job created by the loop at the bottom of this script

	# Read in base ASCII file that is datum WGS1984-UTM, 250 m res, and sized to WTplusBuffer
	coastdist.asc = read.asc('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/coastdistwtplusbuffer250m.asc')

	# Establish row/column positions for 250 m grid
	base.asc = coastdist.asc
	base.pos = as.data.frame(which(is.finite(base.asc), arr.ind = T)) 

	# Add east and north to this data frame of row/column positions
	base.pos$east = getXYcoords(base.asc)$x[base.pos$row]
	base.pos$north = getXYcoords(base.asc)$y[base.pos$col]

	# Convert east-north to lat-long
	### EPSG 32755 WGS1984 Lat Long
	### EGSG 4326 WGS1984 UTM
	tout = as.data.frame(spTransform(SpatialPoints(cbind(base.pos[3:4]),proj4string=CRS("+init=epsg:32755")), CRS("+init=epsg:4326")))
	base.pos$lon = tout[,1]
	base.pos$lat = tout[,2]

	# Write out a .csv file of row/col positions east/north and lat/long for WT plus buffer
	write.csv(base.pos,file=paste(in.dir,'WTplusBuffer.rowcol.eastnorth.latlong.csv',sep=''), row.names=F)

# Establish a list of all files within the uppermost directory, including the full path
	
tfiles = list.files(in.dir,pattern='.grid.gz',recursive=TRUE,full.names=TRUE)

# Establish a temporary directory to write out the .sh files and .Rout files to
# This keeps your workspace more orderly and should the script fuck up all error logs will be stored in one location

tmp.pbs = '/home1/99/jc152199/tmp.pbs/'
dir.create(tmp.pbs)
setwd(tmp.pbs)

# Input location of the script to be run by this batch script

script.file = '/home1/99/jc152199/MicroclimateStatisticalDownscale/AWAPASCII/RawAWAP/ascii.res.5k.to.250m.2run.r'

# Begin a loop to create an .sh  file for each AWAP layer to be resample and then submit the jobs to the HPC

for(tfile in tfiles) 

	{
	
	tname = basename(tfile) #Get the just the name of the ASCII file from the file list
	
	tname = gsub('.grid.gz','.sh', tname) #Use tname to create the name for the .sh file
	
	arg.tfile = paste('tfile="',tfile,'" ',sep='') # Define the argument for tfile 
	
	zz = file(tname,'w') # This command creates a file with the name 'tname' The 'w' command opens the file and writes the following lines of text into it
		
		cat('##################################\n',file=zz) # First blank line of .sh file ### Note the use of '\n', any text following this will be on the next line
		cat('#!/bin/sh\n',file=zz) # First command line for HPC
		cat('cd $PBS_O_WORKDIR\n',file=zz) # Sets the directory to write out summary files as the above directory 'tmp.pbs'
		cat("R CMD BATCH '--args ",arg.tfile,"' ",script.file,' ',gsub('.sh','',tname),'.Rout \n',sep='',file=zz) # This line defines the argument, says where script to run is located, and defines the name for the R.out file
		cat('##################################\n',file=zz) # Last blank line of .sh file
	
	close(zz) # Closes the .sh file that you've just written
	
	# The follow line submits the job to the HPC, don't comment it out until you're ready to run all the jobs
	#system(paste('qsub -m n ',tname,sep=''))

}

# Done